rm(list=ls())
library(foreign)
library(Zelig) #requires Zelig 3.5.5 for compatibility
                #available at https://cran.r-project.org/src/contrib/Archive/Zelig/

#number of imputed datasets
m <- 10 

#delete next line
setwd("C:/Users/Martin Steinwand/Documents/Project Maritime Boundaries/Analysis/R&RJCR/ReplicationPackage")
#uncomment next line, update your path
#setwd(<Your replication directory>)

#load raw data
alldat <- read.dta("Boundary Data_work.8-22-15.dta")
attach(alldat)
if(id2[2]-id2[1]!=0)cat("Data not sorted \n")
if(year[2]-year[1]!=1)cat("Data not sorted \n")
alldat$difshare <- -(alldat$difshare-100)
workdat <- data.frame(subset(alldat, ongoing2==0))
detach(alldat)

#load imputed data
load(file="BDdat8-2-12y5dum.NoOngoing.RData")

#create new variables
newdat <- transform(newdat, ddema=as.numeric(politya>=7))
newdat <- transform(newdat, ddemb=as.numeric(polityb>=7))
newdat <- transform(newdat, autoca=as.numeric(politya<=-5))
newdat <- transform(newdat, autocb=as.numeric(polityb<=-5 & polityb>=-10))
newdat <- transform(newdat, jointdem=as.numeric(ddema==1 & ddemb==1))
newdat <- transform(newdat, jointautoc=as.numeric(autoca==1 & autocb==1))
newdat <- transform(newdat, demautoc=as.numeric((ddema==1 & autocb==1) | (ddemb==1 & autoca==1)))
newdat <- transform(newdat, gdpdif=abs(gdpa2 - gdpb2)/(gdpa2 + gdpb2))
newdat <- transform(newdat, dependsingle=as.numeric((dependa==1 & dependb==0)|(dependb==1 & dependa==0)))
newdat <- transform(newdat, dependboth=as.numeric(dependa==1 & dependb==1))
newdat <- transform(newdat, coastdif=abs(coastla-coastlb)/(coastla+coastlb))
newdat <- transform(newdat, dcap=abs(capa-capb))
newdat <- transform(newdat, lcoastdif = log(as.numeric(coastdif<=0)*.001+as.numeric(coastdif>0)*coastdif))
newdat <- transform(newdat, vetotot=(vetoa + vetob)/2)
newdat <- transform(newdat, vetodif=abs(vetoa - vetob))
newdat <- transform(newdat, tsp1=workdat$tsp1)
newdat <- transform(newdat, tsp2=workdat$tsp2)
newdat <- transform(newdat, tsp3=workdat$tsp3)
newdat <- transform(newdat, totag=workdat$totag)
newdat <- transform(newdat, depend=as.numeric(dependsingle==1 | dependboth==1))
newdat <- transform(newdat, los=as.numeric(lossingle==1 | losboth ==1))
newdat <- transform(newdat, mid3y=workdat$mid3y)
newdat <- transform(newdat, mid5y=workdat$mid5y)
newdat <- transform(newdat, mid10y=workdat$mid10y)
newdat <- transform(newdat, hmid3y=workdat$hmid3y)
newdat <- transform(newdat, hmid5y=workdat$hmid5y)
newdat <- transform(newdat, hmid10y=workdat$hmid10y)
#additional transformations
newdat <- transform(newdat, difshare=workdat$difshare)
newdat <- transform(newdat, difshare.sc=(workdat$difshare*.01*(189-1)+.5)/189)#rescaled for beta model
            #rescaling according to Smithson and Verkuilen (2006)
newdat <- transform(newdat, diffrac=workdat$difshare/100) #for fractional logit
newdat <- transform(newdat, difdisc=as.numeric(workdat$difshare>=50))
newdat <- transform(newdat, difdisc2=as.numeric(workdat$difshare>0))
#create interactions
newdat <- transform(newdat, JdemOildis=jointdem * workdat$oildisc_f)
newdat <- transform(newdat, JdemOilprod=jointdem * workdat$oilprod_f)
newdat <- transform(newdat, JdemTrade=jointdem * workdat$tradetot1000_f)
newdat <- transform(newdat, JdemFish=jointdem * workdat$fishtot_f)
newdat <- transform(newdat, DcapOildis=dcap * workdat$oildisc_f)
newdat <- transform(newdat, DcapOilprod=dcap * workdat$oilprod_f)
newdat <- transform(newdat, DcapTrade=dcap * workdat$tradetot1000_f)
newdat <- transform(newdat, OilProdDisc=workdat$oildisc_f*workdat$oilprod_f)
newdat <- transform(newdat, oilone=workdat$oilone)
newdat <- transform(newdat, oiltwo=workdat$oiltwo)
newdat <- transform(newdat, oileither=workdat$oileither)#selection equation
newdat <- transform(newdat, oilproddisc=workdat$oileither * workdat$oildisc_f)#selection equation
#now add the dummies for directed oil discovery
newdat <- transform(newdat, oildiscboth = as.numeric(workdat$oildisa==1 & 
                    workdat$oildisb==1))
newdat <- transform(newdat, oildiscone = as.numeric((workdat$oildisa==1 |
                    workdat$oildisb==1) & workdat$oildisa!= workdat$oildisb))
newdat <- transform(newdat, oildiscany= oildiscone+oildiscboth)
newdat <- transform(newdat, totdif=workdat$totdif)
newdat <- transform(newdat, ltotdif=log(workdat$totdif+1))

#model 1
sform.strip <- formula(agreement ~ oildisc_f + oileither+ oilproddisc+
            tsp1 + tsp2 + tsp3 +
            y5dum1 + y5dum2+ y5dum3 + y5dum4 + y5dum5 + 
            y5dum6 +y5dum7 + y5dum8 +y5dum9 -1)

s.modstrip <- zelig(sform.strip, 
            data=newdat$imputations, model="probit", cite=FALSE)
summary(s.modstrip)

#model 2
sform <- formula(agreement ~ jointdem + oildisc_f + oileither+ oilproddisc+
            fishtot_f + tradetot1000_f + los + totag+ dependsingle + dependboth + 
             terrdisp +  voteshare + gdpdif + dcap + tsp1 + tsp2 + tsp3 +
            y5dum1 + y5dum2+ y5dum3 + y5dum4 + y5dum5 + y5dum6 +y5dum7 + y5dum8 +y5dum9 -1)

s.mod <- zelig(sform, 
            data=newdat$imputations, model="probit", cite=FALSE)
summary(s.mod)

#model 3
sform <- formula(agreement ~ jointdem + oildiscboth + oildiscone + oileither+ oilproddisc+
            fishtot_f + tradetot1000_f + los + totag+ dependsingle + dependboth + 
             terrdisp +  voteshare + gdpdif + dcap + tsp1 + tsp2 + tsp3 +
            y5dum1 + y5dum2+ y5dum3 + y5dum4 + y5dum5 + y5dum6 +y5dum7 + y5dum8 +y5dum9 -1)


s.mod <- zelig(sform, 
            data=newdat$imputations, model="probit", cite=FALSE)
summary(s.mod)
